% function [MatDiag, DiagExact] = testRho
if(1)
  ham2d;
end


if(1)
  beta = 500;
  % mu = SD(64);
  mu = SD(16);
  RhoExact = real(2d0 * EF * diag( 1./ (exp(beta*(D-mu))+1d0) ) * EF');
end

if(1)
  nLevel = 4;

  [MatDiag, SampleList, RKTree] = ConstructHMatrix( RhoExact, nLevel );
  DiagExact = reshape(diag(RhoExact),N,N);
  err = norm(MatDiag - DiagExact, 'fro') / norm(DiagExact, 'fro')
end

if(0)
  display('Checking error for each level');
  GRD = reshape(1:N*N, N, N);
  ErrBlock = cell(nLevel, 1);
  for iLevel = 1 : nLevel
    display(iLevel);
    ErrBlock{iLevel} = zeros( RKTree{iLevel}.nRKBlock, 1 );
    for iRKBlock = 1 : RKTree{iLevel}.nRKBlock
      Source = RKTree{iLevel}.RKSource(:, :, iRKBlock);
      Target = RKTree{iLevel}.RKTarget(:, :, iRKBlock);
      SourceInd = reshape( GRD(Source(1,1):Source(1,2), ...
	Source(2,1):Source(2,2)), [], 1 );
      TargetInd = reshape( GRD(Target(1,1):Target(1,2), ...
	Target(2,1):Target(2,2)), [], 1 );
      BlockExact = RhoExact( TargetInd, SourceInd ); 
      BlockApprox = RKTree{iLevel}.RKBlockU{iRKBlock} * ...
	(RKTree{iLevel}.RKBlockV{iRKBlock})';
      ErrBlock{iLevel}(iRKBlock) = normest( BlockExact - BlockApprox );
    end
    mean( ErrBlock{iLevel} )
    sqrt( var(ErrBlock{iLevel}) )
  end
end

if(1)
  display('Checking 2-norm for each level');
  GRD = reshape(1:N*N, N, N);
  RhoRKExactAll = zeros(N*N, N*N);
  RhoRKAll = zeros(N*N, N*N);
  Err2Big = zeros( nLevel, 1 );
  for iLevel = 1 : nLevel
    display(iLevel);
    RhoRKExact = zeros(N*N, N*N);
    RhoRK = zeros(N*N, N*N);
    for iRKBlock = 1 : RKTree{iLevel}.nRKBlock
      Source = RKTree{iLevel}.RKSource(:, :, iRKBlock);
      Target = RKTree{iLevel}.RKTarget(:, :, iRKBlock);
      SourceInd = reshape( GRD(Source(1,1):Source(1,2), ...
	Source(2,1):Source(2,2)), [], 1 );
      TargetInd = reshape( GRD(Target(1,1):Target(1,2), ...
	Target(2,1):Target(2,2)), [], 1 );
      RhoRKExact(TargetInd, SourceInd) = ...
	RhoExact( TargetInd, SourceInd );
      RhoRK(TargetInd, SourceInd) = RKTree{iLevel}.RKBlockU{iRKBlock} * ...
	(RKTree{iLevel}.RKBlockV{iRKBlock})';
    end
    Err2Big(iLevel) = normest( RhoRKExact - RhoRK );
    Err2Big(iLevel)
    RhoRKExactAll = RhoRKExactAll + RhoRKExact;
    RhoRKAll = RhoRKAll + RhoRK;
  end
  Err2 = normest(RhoRKExactAll - RhoRKAll)
end
